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Abstract. The observed accelerated cosmic expansion can be a signature of fourth -order 
gravity theories, where the acceleration of the Universe is a consequence of departures from 
Einstein General Relativity, rather than the sign of the existence of a fluid with negative 
pressure. In the fourth -order gravity theories, the gravity Lagrangian is described by an 
analytic function f(R) of the scalar curvature R subject to the demanding conditions that 
no detectable deviations from standard GR is observed on the Solar System scale. Here we 
consider two classes of f(R) theories able to pass Solar System tests and investigate their 
viability on cosmological scales. To this end, we fit the theories to a large dataset including 
the combined Hubble diagram of Type la Supernovae and Gamma Ray Bursts, the Hubble 
parameter H(z) data from passively evolving red galaxies, Baryon Acoustic Oscillations 
extracted from the seventh data release of the Sloan Digital Sky Survey (SDSS) and the 
distance priors from the Wilkinson Microwave Anisotropy Probe seven years (WMAP7) data. 
We find that both classes of f(R) fit very well this large dataset with the present - day values of 
the matter density, Hubble constant and deceleration parameter in agreement with previous 
estimates; however, the strong degeneracy among the f(R) parameters prevents us from 
strongly constraining their values. We also derive the growth factor g = d In 5/d In a, with 
5 = 5pM I ' Pm the matter density perturbation, and show that it can still be well approximated 
by g{z) oc Q,m{zV ' ■ We finally constrain 7 (on some representative scales) and investigate 
its redshift dependence to see whether future data can discriminate between these classes of 
f(R) theories and standard dark energy models. 
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1 Introduction 

To explain the present-time accelerated expansion of the Universe and the temperature 
anisotropy pattern of the Cosmic Microwave Background (CMB) radiation (e.g. [1, 2]), 
the current concordance cosmological model relies on the properties of two so-called "dark 
components:" Dark Matter (DM) and Dark Energy (DE). On the one hand, the need for 
the former - a non-relativistic, weakly-interacting, non-baryonic matter - matches today's 
requests for an extension of the standard model of particle physics: for instance, there are 
several proposals for DM candidates in the framework of supersymmetry [3-5]. However, 
different approaches have also been followed - for example MoND (modified Newtonian dy- 
namics) [6-8] and conformal gravity [9, 10]. On the other hand, DE still represents a difficult 
theoretical challenge (see [11] for a different perspective). Indeed, there are strong fine-tunings 
problems, regardless of whether one interprets DE as a cosmological constant A, whose tiny 
but non-zero value is not supported by any geometrical symmetry, or whether one considers 
DE as the vacuum energy of a quantum field. In this case, the discrepancy between its actual 
measured value and the theoretical estimate is bigger than a hundred orders of magnitude. 

To remove the DE problem, a different approach has recently started to be investigated. 
It is based on the consideration that A was originally proposed as a constant term in the left- 
hand (geometric) side of Einstein's field equations. By generalising this approach, one can 
argue whether it is possible to reproduce the current cosmic accelerated expansion by adding 
a non-constant time-dependent term in Einstein's tensor. Actually, the effort of modifying 
and generalising the Hilbert-Einstein action of gravity dates back to just few years after 
Einstein's seminal papers (e.g. [12] for a historical review), and it has been also proposed 
by Starobinsky [13] in order to explain the cosmic inflation in the early Universe. This idea 
has again been suggested nowadays to correctly describe the current accelerated expansion 
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of the Universe without any exotic fluid [14-19]. These modified-action theories of gravity 
are widely known as f(R) theories, because in the gravity Lagrangian the Ricci scalar R is 
replaced by the addition of a function f(R), which can be, in principle, generic. 

However, GR is a well-tested theoretical framework, at least with respect to the Solar 
System scale of distances. Therefore, any f(R) theory which attempts to solve the late-time 
acceleration problem has to face the Solar System tests of gravity. Recently, two models 
carefully designed to pass the local gravity tests but still providing an accelerated cosmic 
expansion have been proposed [20, 21]. In this work, we analyse and test these viable f(R) 
models. Specifically, we fit these theories against a large dataset which includes the combined 
Hubble diagram of Type la Supernovae and Gamma Ray Bursts, the Hubble parameter H (z) 
data from passively evolving red galaxies, Baryon Acoustic Oscillations extracted from the 
seventh data release of the Sloan Digital Sky Survey (SDSS) and the distance priors from the 
seven years data of the Wilkinson Microwave Anisotropy Probe (WMAP7). Furthermore, we 
derive the growth factor g = din 5 /din a, with 5 = §pm/pm the matter density perturbation. 

2 f(R) theories 

Being a straightforward generalization of Einstein GR, fourth order gravity theories have 
been investigated almost as soon as the original Einstein theory appeared. It is, therefore, 
not surprising that the corresponding field equations and the resulting cosmology have been 
so widely discussed in the literature. We here first review the basic of f(R) theories and then 
introduce the two classes of f(R) models we will investigate in this work. 

2.1 f(R) cosmology 

In the framework of the metric approach, the field equations are obtained by varying the 
gravity action 

'f(R) 



2k 

with respect to the metric components. We obtain 



S = I d^x^p^g 



+ C 



M 



(2.1) 



f'R^u - V^f + {Of - i/J 9uv = kT, w , (2.2) 



where R is the scalar curvature, k = 8ttG (in units with c = 1), Cm is the standard matter 
Lagrangian with Ta V the matter stress - energy tensor, and the prime denotes derivative with 
respect to R. Note that, for f(R) = R — 2A, one obtains the usual Einstein equations with a 
cosmological constant A. In the general case, a further scalar degree of freedom, conveniently 
represented by the scalar function <f> = f — 1, is introduced. The trace of Eq.(2.2) gives the 
evolution for the effective field [22] 

with the potential V related to f(R) by 

^ = i(2/ -/'#)/". (2.4) 
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In a spatially flat homogenous and isotropic universe, the time- time component of Eq.(2.2) 
gives the evolution of the Hubble parameter, H = a/a (with the dot denoting derivative with 
respect to cosmic time t), being: 

U 2 | dQnf) f-Rf K PM 

while the trace equation (2.3) may be recasted as a Klein -Gordon relation for eft: 

ip + m4>+ d ^ = y M ■ (2.6) 

Eqs.(2.5) and (2.6) can be rearranged in such a way that a single equation for the Hubble 
parameter is obtained. To this end, it is first convenient to change variable from the time t 
to the redshift z = 1/a — 1 and define the dimensionless curvature £ = R/m 2 , with 

2 _ K 2 p M (r] = 0) _ ^ ooi E A/r _ (VL M h 2 



m2 = ..^v, „ (8315 Mpc) - a (2 7) 

a convenient curvature scale which depends on the present day values of the matter density 
parameter Qm and the Hubble constant h = i/o/(100 km/s/Mpc). Assuming dust as gravity 
source and introducing E = H(z)/Hq, it is then only a matter of algebra to get : 

p 2M _ ^[(l + ^) 3 + (^-^ 2 /)/6] (2R) 
j — m z {l + z)(d€/az) 



dz 2 



dz 



d\nE 



d\a.(l + z) 



m 2 f" d£ 
1 + z dz 



E 2 



2m 2 f - W 
3(1 + zf 



(2.9) 

By inserting Eq.(2.8) into Eq.(2.9), we get a single second order nonlinear differential equation 
for £(2) that can be solved numerically provided f(R) and the initial conditions are given. 
To this end, we first remember that, because we are using the RW metric, it is 



R = 6(H + 2H 2 ) 
so that we get for the present-day values [23] : 



(2.10) 



Rn 



6H 2 (l-q ) , i? 



6#o 3 (Jo-<7o 



-H 2 d/a) and jerk (j 



(2.11) 



with (f/Oiio) the present-day values of the deceleration (q = 
H~ s 'a /a) parameters. It is then straightforward to get : 

'£(z = 0) = (6/n M )(l-q ) 

< 

k d£/dz(z = 0) = (6/n M )(j -q -2) 

As a final remark, we note that, because of the definition of £, Eq.(2.9) is a single fourth- order 
nonlinear differential equation for the scale factor a(t) so that we need to know the values of 
the derivatives up to the third order to determine the evolution of a(t). This explains why 
the jerk parameter jo is required. 



- 3 - 



2.2 f(R) models 

A key role in fourth - order gravity is clearly played by the functional expression of f(R). 
Although such a choice is in principle arbitrary, there are some fundamental tests that have 
to be passed in order to get a physically viable theory. First, at the Solar System scales, 
f(R) theories introduce a scalar degree of freedom which couples to matter and originates a 
long range fifth force that can lead to incorrect PPN parameters (see, e.g., [24-27] and refs. 
therein). As a possible way out, one can invoke a chamaleon effect [28, 29] and tailor the 
f(R) expression in such a way to give rise to a mass squared term which is large and positive 
in high curvature environments [30-33]. On the other hand, in the early Universe, one wants 
to recover the standard GR in order to preserve the agreement with the nucleosynthesis 
constraints. Among the possible choices left out by the above conditions, we will consider 
here two popular classes of f(R) models which we briefly describe in the following. 

After working out the above conditions for the viability of a f{R) theory, Hu & Sawicki 
first proposed the following functional expression [20] : 

f(R) = R- m \ Cl{R/ ™X (2.12) 
Jy ' l + c 2 (R/m 2 ) n y ' 

where m 2 is given by (2.7), and (n,c\,C2) are positive dimensionless constants. Note that, 
since f(R = 0) = 0, there is formally no cosmological constant term. Nevertheless, since 

v ftm Cl 2 , c i 2 ( ™ 2 
hm t{R) — m H — — 

we still recover an effective A term in high curvature (m 2 /R — > 0) environments. In particular, 
in the limit Rq » m 2 and Rq » -R*, the expansion rate H, in the early universe, will be 
the same as in ACDM with an effective matter density parameter £lM,eff = 602/ (ci + 6C2) 
that guarantees that the nucleosynthesis constraints are satisfied. 

The HS model is determined by the three parameters (n, 01,02): two of them can be 
fixed in terms of observable quantities. To this end, we first note that evaluating Eq.(2.8) 
at z = gives a relation between the derivatives of / and the present day value of d^/dz. 
Second, we expect that f'(R = Rq) = fm only mildly departs from the GR value fjio = 1 
in order to have an effective gravitational constant G e fj = G/f as close as possible to the 
local one today. We therefore have : 

»Af[l + (Co/iffl ~ m 2 f )/6] = i 

fm — m 2 fo£po ^2 13) 

fm = 1 - e 

with £0 = £(z = 0), fo = f(R = Ro), £ p o = d£/dz(z = 0), and e an additional free parameter. 
Inserting the corresponding expressions for the HS model, Eqs.(2.13) can then be solved to 
express (01,02) as a function of (fijif? qo,jo, e) thus simplifying the choice of the parameters 
when fitting the model to the data. 

Another possibility to satisfy all the constraints hinted at above is offered by the 
Starobinsky proposal [21] : 



f(R) = R + \R* 



i R 2 \ 

1 + 



(2.14) 
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with R* a scaling curvature parameter and (A, n) two positive constants. We will refer to 
this class of f(R) theories as the St model. Note that, even in this case, /(0) = so that 
no actual cosmological constant is present. Nevertheless, an effective one is recovered in the 
high curvature regime as can be seen from f(R 3> R+) ~ R — 2Aoo with Aoo = XR+/2. The 
St model parameters are (n,i?*,A), but it is again convenient to reparametrize the model 
differently. We again resort to Eqs.(2.13) inserting the corresponding expressions of / and /' 
for the St model and solving them with respect to (A, it!*). Note that we thus use the same 
set of parameters to describe both the HS and the St models. 

These two classes of f(R) theories are actually quite similar at both very low and 
very high redshifts since they are both built up by imposing the same constraints on f(R). 
Moreover, they both aim at mimicking the successful ACDM scenario in the late and early 
Universe. In other words, the HS and St models both reduces to the GR + A case in the limits 
of very high and very low curvatures. What makes them different is the way the two extreme 
cases are connected, i.e. how the Universe evolves from the present day phase dominated by 
A to the early matter epoch. 

For completeness, we finally remind the reader that the two models we are considering 
here are not the only viable ones; other possible examples are given in [35, 36]. It is, moreover, 
possible to design f(R) models that can also provide an inflationary expansion in the very 
early Universe [37, 38]. However, all these other cases share many similarities with the HS 
and St models so that we are confident that exploring these two classes of fourth -order 
gravity theories should provide some general conclusions on their viability. 

2.3 The scale factor evolution 

To obtain the numerical solutions of Eqs.(2.8) and (2.9) more efficiently, we resort to an ana- 
lytical approximation. We note that, for both the HS and St models, the gravity Lagrangian 
is excellently approximated by the GR + A one by construction. Therefore, unless we are 
interested in the very early universe (where both f{R) models have a vanishing A term), we 
can expect that the dimensionless Hubble parameter, E(z) = H(z)/Hq, is close to ACDM 



E( z » ZA ) ~ e x (z) = yjn M;eff (i + z) 3 + n 7 (i + zY + (i - n M>eff - n r ) (2.15) 

where z\ is a characteristic redshift marking the transition to the ACDM regime, 7 is the 
radiation density parameter and 



n 



M,eff 



6c 2 



ci + 6c 2 



l ~^ fOT St 



for HS 

(2.16) 



is the effective matter density parameter as inferred from the effective A term introduced in 
this limit by the HS and St models. 

On the other hand, from the point of view of the background evolution, every modi- 
fied gravity theory is equivalent to a dark energy (DE) model with a suitably reconstructed 
equation of state. Moreover, for a large class of DE theories, the equation of state is well ap- 
proximated over a large redshift range by the phenomenological Chevallier - Polarski - Linder 
(CPL) [39, 40] parametrization, w(z) = wq + w a (l — a), leading to 
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e 2 cpl {z) = n M (i + zf + (i - n M )(i + z) 3(i+-o+^) exp 

for the dimensionless Hubble parameter so that, for (wo,w a ) = (—1,0), E A (z) is recovered. 

Motivated by these considerations, for both the HS and St models, we have therefore 
fitted the numerical solution of Eqs.(2.8) and (2.9) with the following ansatz : 

( £(z)E C PL(z,n M ) + [1 - £{z)]E A (z,n M ) z < z A 
E{z) = I (2.18) 

[ E A (z,flM,eff) Z>Z A 

where 

3 

£(z) = Y,e l (z-z A ) i (2.19) 
i=i 

is an interpolating function with and z A fitting parameters. This approximating function 
excellently reproduces the numerical solution whatever the model parameters (Qm, Qo,jo, n , e) 
are for both the HS and St models with a rms error which is far lower than 0.1% over the full 
redshift range (0, 1000). A simple look at the f(R) functions makes quite easy to understand 
why this happens. Indeed, both f(R) functions converges to f(R) = R — 2A e ff as z increases 
so that it is not surprising that, after a critical value z A when f(R) is indistinguishable from 
the ACDM Lagrangian, the Hubble parameter becomes exactly the same with the value of 
the effective cosmological constant depending on the model parameters. On the other hand, 
at very low redshift, the CPL parametrization does an excellent job in mimicking E(z « 1) 
as expected. The interpolating function £ (z) then only smooths the transition between the 
two regimes in an efficient way. 

A subtle remark is in order here concerning the value of Om- Indeed, while for z < z A , 
we use the actual matter density parameter (and neglect radiation for simplicity) , the effective 
one enters the z > z A approximation. Therefore, a discontinuity in z A is formally present in 
our approximation. Actually, it is easy to show that, for all reasonable model parameters, Qm 
and VtM,eff are almost perfectly equal so that the discontinuity can not be detected at all and 
E(z) is, to all extents, a continuous function. On the other hand, it is worth stressing that 
(wo,w a , ei, e2, e3, z A ) depend in a complicated way on the f(R) model parameters. Actually, 
by numerical attempts, we discovered that, within a very good approximation, (ei,e2,es) 
are the same for all the models considered so that the four remaining quantities (qo,jo,n,e) 
collapse into three parameters only, namely (wo,w a ,z A ). We can therefore anticipate that 
strong degeneracies among them will take place since different sets will give rise to the same 
E(z). Moreover, z A increases with e because the smaller is e, the closer is Jrq (and hence 
/') to 1, i.e. the smaller e is the quicker is the convergence to the ACDM Lagrangian. As 
a consequence, the smaller is e, the narrower will be the redshift range where E{z) departs 
from the concordance model. 

3 Constraining f(R) models 

The HS and ST models provide an accelerated expansion in a matter only universe because, 
by construction, they mimic the successful ACDM scenario in both the low and high redshift 
regimes. Nevertheless, fitting the model to the data is still of significant help to constrain 
their wide parameter space and pin down their predictions on other not fitted quantities. In 



3w a z 
l + z 



(2.17) 
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order to lift the degeneracies among the model parameters, one can not rely on low redshift 
probes only, but has to add further data tracing higher z or based on quantities related in 
a different way to the Hubble parameter. Such considerations motivated our choice of the 
dataset described below. 

3.1 Observational data 

We use the 557 SNela in the Union2 [1] sample to probe the evolution of the low redshift 
universe (over the range 0.015 < z < 1.37). Neglecting systematics, the SNela likelihood 
term will simply read 

C S NeIa(p) OC exp [-XSNelaiP) l 2 \ C 3 - 1 ) 

with 

NsNela r I \ / \ 1 2 

Vobs(Zi) ~ fi th {Zi,p) 



XSNelaiP) = 



(3.2) 

i=l 

where Oi is the error on the observed distance modulus fj, b s {zi) for the i-th object at redshift 
Zi and the theoretical distance modulus is given by 



lkh(z,p) = 25 + 5 log 
with r(z) the dimensionless comoving distance 



Q 

— (l + z)r(z,p) 



(3.3) 



/ \ r dz' , 

'• (Z ' P) = X WT)- (3 ' 4) 

and p the set of model parameters for the assumed f(R) class of theories. 

The HS and St models differ from each other in the intermediate redshift regime and we 
need to trace the Hubble diagram in the matter dominated era to discriminate between them. 
Thanks to the enormous energy release that makes them visibile up to z ~ 8.2, Gamma Ray 
Bursts (GRBs) stand out as good candidates to this task [41, 42]. We therefore use the GRBs 
Hubble diagram as recently derived in [43] (see also [44]) for the catalog of 115 GRBs in [45]. 
In order to use a model independent calibration of the GRBs scaling relations (but see [46] for 
the limits of these assumption), we use the results in [43] obtained using the local regression 
method (see [42] and refs. therein for details). From the sample, we remove the objects with 
z < 1.4 in order to avoid correlations, difficult to quantify, with the SNela data introduced 
by the calibration method. As a result, we end up with 64 objects with 1.48 < z < 5.60, so 
probing the Hubble diagram deep into the matter dominated era. 

Note that, when using GRBs, the likelihood term will be the same as Eq.(3.1), but 
with the denominator in Eq.(3.2) changed to (cr? + of nt ) 1//2 with <Ji n t the unknown intrinsic 
scatter of the GRBs around the theoretical Hubble diagram. This term may come out from 
different sources. First, the empirical GRBs Hubble diagram has been obtained by averaging 
the distance modulus of each object as inferred from multiple scaling relations so that, if 
an object is an outlier for one of the correlations used, its distance modulus will be shifted 
from the true one. Moreover, the calibration procedure may also introduce some bias due 
to neglecting any redshift evolution of the scaling relations coefficients which can not be 
excluded given the wide redshift range probed. Therefore, we leave Oi nt as an unknown 
parameter 1 and marginalize over it in the likelihood analysis. 



1 Note that a similar term is also present for SNela, but it is estimated to be ai nt = 0.15 and yet included 
into the error <t; provided in the Union2 dataset. 
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While both SNela and GRBs are based on the concept of standard candles, an alterna- 
tive way to probe the background evolution of the universe relies on the the use of standard 
rulers. A nowadays widely used example is represented by the Baryonic Acoustic Oscillations 
(hereafter, BAOs) which are related to the imprint of the primordial acoustic waves on the 
galaxy power spectrum. In order to use BAOs as constraints, we follow [47] by first defining : 



r s (zd) 
D V (z) 



(3.5) 



with Zd the drag redshift (computed using the approximated formula in [48]), r s (z) the 
comoving sound horizon given by : 



r s(z) 



da 



VSJo a 2 H(a)y/l + {3/4)n b /n 7 

and Dy{z) the volume distance defined by [49] : 



(3.6) 



D v (z) 



cz 
H[z) 



Dl(z) 
l + z 



1/3 



(3.7) 



We finally constrain the model parameters by introducing the likelihood function 



£bao(p) 



exp [-(Dbao ■ C 



BAO 



D 



T 

BAO 



)/2] 



9-7rlP~ 1 1 1/2 



(3- 



with T> T BAO 



(^a! ~~ ^a2>^a35 ~~ ^035) an d Cbao the BAO covariance matrix. The values 
of d z at z = 0.20 and z = 0.35 have been estimated by [47]; we also use the values of the 
observed d z and C bao provided by [47] . 

Both the Hubble diagram and the BAO are distance related quantities so that they probe 
the integrated expansion rate. Thus, the details of H{z) are partially smoothed out so that 
severe degeneracies arise. In order to alleviate this problem, one should rely on direct esti- 
mates of H{z) which can be obtained by noting that, from the relation dt/dz = —(l + z)H(z), 
a measurement of dt/dz at different z gives, in principle, the Hubble parameter. This differen- 
tial age method [50] works best using fair samples of passively evolving galaxies with similar 
metallicity and low star formation rate so that they can be taken to be the oldest objects at 
a given z. Stern et al. [51] have then used red envelope galaxies as cosmic chronometers de- 
termining their ages from high quality Keck spectra and applied the differential age method 
to estimate H(z) over the redshift range 0.10 < z < 1.75 [52]. We use their data as input to 
the following likelihood function : 



£h(p) 



with 



Xh 



£ 

i=l 



exp [-x 2 H (p)/2] 
(2vr)^/ 2 |C^ 1 | 1 / 2 ' 



Hobs(zi) - H(zi,p) 



(3.9) 



(3.10) 



where the sum is over the Mr = 11 sample points and Ch is the diagonal covariance matrix. 
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The set of data described so far probes the evolution of the universe only during the 
late (dark energy dominated) era and the matter epoch, but it tells us nothing about the 
early Universe. This problem can be healed by adding the CMBR data to have a picture 
of the Universe at recombination [z ~ 1100). Komatsu et al. [53] have shown that most of 
the information in the WMAP power spectrum may be summarized in the so called distance 
priors, i.e. constraints on : (1) the redshift zls °f the last scattering surface (computed with 
the approximated formula in [54]), (2) the acoustic scale [55-57] : 

ir(c/H )r(z LS ) , . 

i A = - f x , (3.11) 

r s {zLS) 

and (3) the shift parameter [55-57] : 

K = ^/n^r(zLs) . (3.12) 
We can then define the likelihood fucntion 

r exp [-(Dca/b • Ccm B ■ gcMgVg] , . 

Wp) = W\cbU 1 » ' ( ] 

with ~Dqmb the vector with the values of the differences between the observed and the 
theoretically predicted distance priors and Ccmb the corresponding covariance matrix. We 
rely on the seventh data release of the Wilkinson Microwave Anisotropy Probe (WMAP7) 
[2] to set the observed distance priors and their covariance matrix. 

3.2 Likelihood analysis 

In the Bayesian approach to model testing, the parameter space of a given model is con- 
strained by examining the region where a user -defined likelihood function £(p) takes non 
negligible values. In particular, the best fit parameters are the set Pbf maximimizing £(p), 
while the constraints on the i-th quantity pi are obtained by marginalizing C over all the 
parameters but pi itself. Mathematically, one defines : 

£i(Pi) = J dpi... J dp^i j dp i+i ...J dp n £(p) (3.14) 

and estimates the median value and the 68 and 95% confidence ranges for pi from Ci(pi). 
As a general remark, we warn the reader that, in the context of Bayesian statistics, the best 
fit model represents the most plausible model in an Occam's razor sense given the data at 
hand. However, in a Bayesian context, the best fit parameters individually do not necessarily 
have to be probable, but rather they must have a high joint probability density that might 
occupy only a small volume of the parameter space. This situation can arise if the best fit 
solution does not lie in the bulk of the posterior probability distribution. Such a case may 
often occur when the posterior is non-symmetric in a high dimensional space so that the 
volume can dramatically increase with the distance from the best fit solution. In this case, 
the best fit solution for each parameter could easily lie outside the bulk of the individual 
posterior distribution for pi obtained by marginalizing over the other parameters. This is 
indeed what happens for our models so that we have preferred to remind the reader that 
this somewhat counterintuitive outcome is not a problem, but rather a common feature in 
statistics in multidimensional spaces. 
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The likelihood £(p) takes the available data and the prior information on the model 
into account. In order not to bias our analysis, we include only flat priors on (n, loge) as 
explained later. We can therefore set : 



£(p) 



exp [-(h s -h) 2 /2a 2 s ] 



(3.15) 



2ira 2 



where the first factor accounts for the recent determination of the Hubble constant from local 
distance calibrators by the SHOES collaboration which give (hs,(Js) = (0.742,0.036) [58], 
while (j) stands for SNeIa,GRB, H, BAO,CMB according to the definitions given in the 
previous subsections. 

As a preliminary task, one should ask what the number of parameters to constrain 
is. First, we note that, in order to solve Eqs.(2.8) and (2.9), we must know the values of 
five quantities, namely (&Mi Qo, jo, n, log e), where we have moved to a logarithmic value 
for e since s is expected to be very small. A sixth parameter is the Hubble constant h 
which enters through both the prior in Eq.(3.15) and the analysis of the H{z) data. When 
including GRBs in the total likelihood, we need to estimate also the intrinsic scatter ai n t in 
order to get the errors on the individual GRBs distance moduli. A further quantity is the 
physical baryon density ujf, = Qt,h 2 entering directly through the WMAP7 distance priors and 
indirectly through its effect on the determination of (zd,ZLs)- Finally, the radiation term 
can not be neglected in the early universe so that we add this contribution when using the 
distance priors as constraints. Summarizing, the full set of parameters we consider reads : 



while we set tot = 2.267 x 10~ 2 [53]. We use flat priors on all the parameters choosing large 
enough ranges to avoid biasing the likelihood analysis by any theoretical prejudice. We have, 
however, make two exceptions to this rule forcing n > 1 so that the chamaleon effect can 
take place [33] and only considering models with : 



where the lower limit is simply set because of loss of sensitivity for smaller values, while the 
upper one is obtained by considering that, in order not to have significant deviations from 
the Newtonian gravitational potential on galactic scales, the constraint [20] 



with v max the maximum rotation velocity of the stars, must be satisfied. This would give us 
loge < —6, but to be conservative we will extend the above range to include higher velocity 
systems (such as clusters of galaxies). Moreover, since zls is actually not very sensitive 
to the model parameters (so that zls — 1090 could be used for all the models within an 
excellent approximation), we do not expect to put strong constraints on J7 7 and we can 
marginalize over it. Finally, we also marginalize over Gi n t since this nuisance quantity has 
not a straightforward physical interpretation being most related to the calibration of GRBs 
scaling relations rather to any physical effect. 

In order to explore the eight dimensional parameter space, we use a Monte Carlo Markov 
Chain (MCMC) algorithm running three parallel chains of ~ 25000 points each and checking 



p= (ttM,Q*f,h,qo,jo,n,loge,a in t) 



13 < loge < -3.0 




(3.16) 
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convergence with the usual Gelman - Rubin statistics [59]. After cutting out the burn in 
phase and thinning to avoid spurious correlations, the final merged chain also allows us to 
evaluate the constraints on some derived quantities of interest. To this end, we evaluate the 
parameter density distribution 0(p) along the chain and study the histogram of the values 
thus obtained to infer the corresponding median and confidence ranges. 

3.3 Results 

Since we have used the same parameters for defining both the HS and St models, it is worth 
discussing the constraints on these quantities in parallel thus also checking if they depend on 
the assumed f(R) functional expression. In particular, the best fit values turn out to be : 

(n M ,h,q ,j ,n,loge) = (0.276,0.719,-0.585,0.995,1.53,-5.95) , 
for the HS model, and 

(n M ,h,q ,jo,n,loge) = (0.273,0.725,-0.599,-0.058,1.34,-10.02) , 

for the St one. How well the best fit models reproduce the data may be quantitatively judged 
noting that, for the HS model, it is (with x 2 = X 2 /d-o.f) : 

x|7Ve/a = 0.99 , X 2 GRB = 1.19 , j& = 2.58 , 
(do.2, 4.35) = (0.1893,0.1137) , (£ A ,K,z LS ) = (302.29, 1.725, 1091.5) , 
while the St best fit model gives : 

X|7v e /a = 0.99 , XGRB = m . Xh = 2-56 , 

(do.2, do.35) = (0.1901,0.1142) , (£ A ,K,z LS ) = (302.36,1.723,1091.5) . 

Comparing with the observed values, we can therefore safely conclude that the both the HS 
and St models are in good agreement with the data. A cautionary note is in order here 
concerning the x 2 values reported above and hereafter. The MCMC code maximizes the 
likelihood (3.15) which is strictly not the same as minimizing the single reduced x 2 entering 
its definition. Moreover, from a statistical point of view, the significance level of the above 
reduced x 2 values can not be estimated from the usual tables since these standard results 
do not take into account systematic errors or intrinsic scatter. It is also worth stressing 
that we have defined above the reduced x 2 computing the number of degrees of freedom 
as Md — Mp with Md (Np) the number of datapoints (parameters). While this is formally 
correct, some caution is needed when high values are obtained. Such large numbers may 
indicate that some of the parameters are actually unnecessary to fit a particular subset of 
the data. For instance, ai n t does not enter at all in the fit to the SNela Hubble diagram and 
the H(z) measurement, but it is nevertheless included in J\f p when normalizing the x 2 for 
these datasets 2 . In particular, the high Xh value is simply due to the inclusion of both Oi nt 
and fi 7 which do not affect at all the fit to this subset. Indeed, the best fit models closely 
follows the H(z) data even if xh ^ s quite large. As a general remark, we therefore warn the 
reader to not overrate the reduced x 2 values. 



2 Note that, although we finally marginalize over both f2 7 and (Tint, they are nevertheless varied in the 
fitting procedure so that they must be included in the Af p value used to normalize the y 2 . 
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-10.03 


-7.57 


H c;q+4-09 +4.91 
°- oy -1.17 -1.81 



Table 1. Constraints on the HS and St model parameters. For each parameter, we give the best fit, 
mean and median values and the 68 and 95% confidence ranges. 

The comparison of the best fit models shows that both models predict the same matter 
abundance and quite similar values for the present day Hubble constant and deceleration 
parameter, while strikingly different results are obtained for the jerk parameter with the HS 
model giving almost the same value as predicted for the concordance ACDM scenario (i.e., 
jo = 1). However, one should better rely on the median value for the jo parameter since 
this latter takes fully into account the shape of the jo distribution. As Table 1 shows, the 
median jo values for the HS and St models are fully with the 68% confidence ranges well 
overlapped. Moreover, the constraints in Table 1 are in good agreement with those in the 
literature (see, e.g., [22, 60-62] and refs. therein). To this regard, it is worth stressing that 
our analysis differs from the previous ones in two important aspects. First, we use a more 
recent compilation of data (both SNela and BAO observational constraints) and add the 
GRB Hubble diagram to probe the matter dominated epoch. Second, we leave all the model 
parameters free to vary thus exploring the parameter space without imposing any limitation 
a priori such as, e.g., fixing Qm (as in [60, 62]). Needless to say, such an approach does not 
come for free, the cost to pay being the rise of severe degeneracies in the 6D parameter space. 
An example is, indeed, provided by the jo distribution for the St model with a best fit value 
strongly different from the median one. However, the ACDM value jo = 1 is well within the 
68% confidence ranges for both the HS and St models so that, from the point of view of the 
cosmographic parameters, the constraints are still too weak to discriminate between the two 
models. Unexpectedly, the constraints on (qo,jo) f° r the St model are much weaker than for 
the HS case. Quite weak constraints are obtained in both cases for n and log e. In particular, 
the confidence ranges for n are strongly asymmetric due to the hard prior n > 1 we have set 
to allow for a chamaleon - like effect to take place. In order to explore why this effect takes 
place, we can first consider the constraints on the original quantities entering the two f(R) 
expressions. For the HS model, we find (median value and 68 and 95% confidence range) : 

loff n - 3 47+ 2 - 01 + 3 - 26 loe co - 2 28+ 2 ' 01 + 3 ' 30 

lOgt-l — -0.91 -1.92 ' 10 & 62 — Z ' Z5 -0.90 -1.91 ' 
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so that Eq.(2.12) effectively reduces to f(R) — R — m 2 c\/c2- The HS model is therefore 
quite similar to ACDM; this explains why the constraints on jo give values narrowly centred 
around the ACDM jo. For the St model, we get : 

logA = 1.50tS;jg±?$ , log(iV^) = -l-73t?:Pil;S . 

so that, in the regime R ~ Ro, the effective Lagrangian reduces to f(R ~ i?o) — -^o + 
XRqk(1 — K 2n ) with k = R+/Rq. Since k << 1, we still have an effective cosmological 
constant term; however, for large values of n and not negligible values of k, the second term 
in parentheses, K 2n , might not be negligible. Therefore, values of jo 7^ 1 are possible and this 
explains the large confidence ranges. 

Finally, as an alternative way to describe the background evolution for the HS and St 
models, we can work out the effective DE equation of state as : 



l + w eff (z) 



2 dhxE(z) VL M 0- + zf 



1 



»A/(1 + Z 

E 2 (z) 



1 3 



3ciln(l + z) E 2 (z) 

The effective EoS for both f(R) models may be anticipated considering that the dimensionless 
Hubble parameter E(z) may be well approximated by Eq.(2.18). Whichever are the model 
parameters, the f(R) function reduces to /a at high z so that the EoS reduces to w e ff(z >> 
1) = —1 identically. On the other hand, at low z, both HS and St models approximate 
well the ACDM evolutionary history so that we expect only small deviations. That this is 
indeed the case is confirmed by the values of wq = w e ff(z = 0) and w\ = dw e ff/dz(z = 0). 
Evaluating these quantities along the chain, we get for the HS model : 

(1 + Wo) X 1(J — — l.y_ 5 5 _ 19 7 , MIlXlU — 2 _ 6 . 2 i 

so that we actually find only tiny deviations from the ACDM values (wq,wi) = (—1,0). On 
the contrary, more freedom is allowed by the St model. In this case, we obtain : 

(^ _i_,„ \ v m- 1 — nm+ - 02 + - 04 „, — n io+3-0 +10.7 
(1 + woJxlU — -U.U1_ 1 6 _ 48 , w\ — — U.12_ 41 _ 88 , 

As it is apparent, an effective cosmological constant is still well consistent with these con- 
straints, but also not negligible deviations from ^o = —1 and a varying EoS are allowed. 
Since the jerk parameter depends on both wq and w\ (see, e.g., [23] for its expression for 
the CPL model), we expect to be able to set only weak constraints on jo when fitting the St 
model to the data, while the opposite conclusion may be drawn for the HS case. 



4 Time related observable quantities 

To constrain the cosmological parameters, one can, in principle, use time -related quantities 
such as the lookback time [63] or the age of old high redshift galaxies (OHRG) [64-66]. 
Unfortunately, at the moment, such a test can give only weak constraints since the data 
are too sparse and affected by large errors and possibly by systematic biases still to be 
investigated in detail. We have therefore not used this test in our analysis, but we check here 
a posteriori whether the HS and St f{R) models are consistent with the age data. To this 
end, we follow [66] and consider the age of the Universe at redshift z : 

f z dz' 

'(*■>>- *y wrmAd (4 ' 1) 
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Figure 1. Estimated incubation time for the HS (left) and St (right) f(R) models for the galaxies in 
the sample given by [66]. Short and long dashed lines refer to the central value and 1 and 2a ranges 
for ti nc as estimated from stellar population synthesis models. 



with tH = 9.78/i Gyr the Hubble time. Given a galaxy at redshift zq, we can get an 
estimate of t(zc) as the sum of the age to of the galaxy (inferred by fitting its observed 
spectrum to stellar population synthesis models) and the incubation time ti nc . This latter 
gives us an estimate of the amount of time between the beginning of structure formation and 
the actual formation time of the galaxy. We then use the OHRG sample assembled in [66] 
from the literature then available and, for each galaxy, estimate the incubation time : 

tinc = tmod(zG) ~ t bs(zG) 

with the error naively estimated as : 



Vine — Y 0- mod + (T obs , 

where the quantities with the underscript mod are estimated evaluating t(zc) along the 
chains. Assuming that the formation redshift is approximately the same for all the galaxies 
in the sample (which is a reasonable approximation) , the incubation time estimated from our 
model should be the same for all the objects in the sample. Moreover, its value should be in 
agreement with t{ nc = 0.8 ± 0.4 Gyr estimated based on stellar population evolutionary time 
[67]. To this end, we therefore compute ti nc from the chains obtained before fitting the HS 
and St models to the full dataset and plot the results in Fig. 1. As it is apparent, the errors 
are still quite large so that a definitive conclusion can not be safely drawn. Nevertheless, 
we note that, for all galaxies but that at the lowest redshift (which is however consistent 
within the 2a error bar), tj nc is in agreement with the estimated value for both the HS and 
St models. Thus, we do not find any age problem for these f(R) theories. 

The age data are difficult to extend to higher redshift because of the prohibitively long 
exposure time needed to get a sufficiently high quality spectra of galaxies at z > 2. The situ- 
ation is, however, better for quasars as it is well illustrated by the case of APM 08729 + 5225. 
With a redshift z = 3.91, this object has an estimated age of 2 - 3 Gyr [68] with a best fit age 
tAPM = 2.1 Gyr and a lower la limit tAPM > 1-7 Gyr [69]. Note that the incubation time is 
here not well defined, but it is likely to be very small. Asking that the age of the Universe at 
z = 3.91 is larger than the age of APM 08729 + 5225 allows to set constraints on dark energy 
models [70-73] and shows that, for many models (including the concordance ACDM), the 
condition t{z = 3.91) = £3.91 > t^pM is n °t easy to fulfill. From our chains, we get : 
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for the HS and St models, respectively. For both cases, the 95% upper limit on £3.91 is 
smaller than the la lower limit of tAPM so that we get the same age problem as for the 
ACDM scenario. Although one can repeat our likelihood analysis explicitly including a prior 
on £3.91, the agreement with the extensive dataset considered above makes us confident that 
the problem could reside in the uncertainties in dating high z quasars. This can lead to 
bias high the estimate of tAPM rather than in a failure of the f(R) theories. As a further 
evidence in favor of this interpretation, one may note that it is actually quite difficult to get 
*3.9i > tAPM for almost all the models in the literature. 

We then consider the present day age of the Universe, to = t( z = 0), obtaining (in Gyr) : 



for the HS and St cases. Both these constraints agree well with previous ones in the literature 
(see, e.g., [2] and refs. therein), but are smaller than the estimated age of old globular clusters 
in the Milky Way, tec — 14 Gyr [74]. This discrepancy is very small and, in fact, a different 
estimate of the globular cluster age, tec = 12.6l 2 ' 6 Gyr [75], removes the discrepancy for 
both f(R) models. 

As a final general remark on these time based tests, we note that the constraints on 
both £3.91 and to turn out to be essentially the same for both f(R) models considered. This 
is a further consequence of how very similar both models actually are. Indeed, for z > 2, 
both f(R) functions are actually well approximated by R — 2A e ff, while the constraints on 
each model parameters simply translate into the same effective A term. As a consequence, 
discriminating among them is quite difficult both with distance and time related quantities 
because the smoothing of the evolution rate due to the integration further cancels out the 
subtle differences in the two f(R) functions. 

5 The growth of perturbations 

As extensively shown here, both the HS and St models provide an evolutionary history that 
can be hardly distinguished from ACDM. This is not surprising given that f(R) reduces to 
/a(.R) for R >> R s with R s a characteristic curvature value depending on the model. Since 
the likelihood analysis points towards model with small values of R s , the condition R > > R s 
is soon fulfilled for most of the redshift range probed (from z = to the last scattering surface 
zls)- In order to discriminate between the HS and St models (and, more generally, between 
f(R) theories and dark energy), one must resort to the observables probing the growth of 
perturbations such as cosmic shear [76-79]. As a preliminary analysis, we will consider here 
some theoretical constraints on quantities related to the growth of structures in the two f(R) 
models we are considering. 

5.1 The growth index 

We start by considering here the growth index of perturbations. In the subhorizon limit, a 
primordial matter perturbation 5 = Spm/pm increases according to [80] : 




Gyr , to = 13.34^;!^ Gyr 



5 + 2H5- 4vrg eJ/ (a, k)p M S = 



(5.1) 



with k the wavenumber and 
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Figure 2. Growth index 7 as function of logfc for the HS (left) and St (right) models. Here and in 
the following plots, the bars refer to the 68% CL with the point denoting the median value. 



G e ff(a, k) 



G l + A(k 2 /a 2 )[f"(R)/f'(R)} 



(5.2) 



f'(R)l + 3(k 2 /a 2 )[f"(R)/f'(R)} 

the scale dependent effective gravitational constant. Note that, differently from standard 
GR, the growth equation now depends on the scale k so that we will have a scale - dependent 
growth of perturbations. Rather than solving Eq.(5.1), we can solve for the growth rate 
g = din 5 /din a. Using the redshift z as variable, Eq.(5.1) leads to: 



dgjz) 
dz 



+ 



1 dlnE 2 (z) 
2dln(l + z) 



(2 + 5) 



g(z) 3n M (l + z) 
1 + z 2 



Geff(k, 



E 2 (z) 



G 







(5.3) 



with the initial condition g{z ~ 100) = 1, i.e., we ask that at very high z the growth rate is 
the same as the one for a matter only universe in GR. If GR holds, a useful parametrization 
for the growth rate is given by [81-83] : 



n AI (i + z ) 

E 2 (z) 



(5.4) 



where 7 is referred to as the growth index. For dark energy models with constant EoS, it is 
7 = 3(w - l)/(6w - 5) which give 7 = 6/11 ~ 0.545 for ACDM. 

For f(R) theories, Eq.(5.3) is scale dependent so that we must check whether the ap- 
proximating formula (5.4) holds at each wavenumber k. As a first result, we have fitted the 
approximated relation : 



log g(z,k) = log 50 + 7!og 



E 2 {z) 



(5.5) 



to the solution of Eq.(5.3) for both the HS and St models and different values of the model 
parameters and wavenumber. We indeed find that Eq.(5.5) fits the numerical solution with 
a rms error which, although increasing with k (up to k ~ 0.3 Mpc -1 ), is always smaller than 
2% (1%) for the HS (St) model. Moreover, go is always unity within less than 0.1% consistent 
with the prediction of Eq.(5.4). That such an approximation for the growth factor g(z,k) 
indeed works is actually not a surprise considering the results in [84-86] which have already 
investigated this issue for the St - like model. We nevertheless note that their approach is 
slightly different since, in [84], the authors integrate an equation for g(z,k) as a function of 
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^m(z), while an equation for 7(2) is derived and solved in [85, 86]. On the contrary, we here 
define 7(2, k) as a fitting parameter to a numerically derived (7(2) vs Qm(z) relation for any 
given k. As a consequence, the different results can not be quantitatively compared. 

It is worth stressing that 7 depends on k so that it is interesting to see whether this 
signature may be used to discriminate between the two classes of f(R) theories and between 
them and dark energy. To take into account errors on the model parameters, we have 
therefore fitted (5.5) to the numerical solution of (5.3) for a given k and for all the points in 
the chain. The results are shown in Fig. 2 for both the HS and St models. Here (and in the 
following plots), we report the 68% confidence range obtained by evaluating 7(logfc) along 
the final chain with the point denoting the median value. Note that the resulting error bars 
are typically asymmetric as a consequence of the non Gaussian distribution of the 7 values. 
Moreover, for the St model, the 7 distribution is actually quite narrow, but some outliers 
generate long tails towards small 7 thus leading to strongly asymmetric error bars. 

There are two remarkable lessons which can be drawn from these plots. First, 7 is well 
determined for log A; < —1.5 for both f(R) theories, while a stronger sensitivity to the model 
parameters takes place for larger k values. Second, the scale dependence of 7 is negligible at 
all for log A; < —1.5 with both models having 7 = 0.55 within an excellent accuracy. On the 
contrary, on smaller scales (i.e., larger k), the St model has still a growth index essentially 
constant, while a significative scale dependence is present for the HS model with 7 becoming 
smaller than the ACDM value. This characteristic features suggest that a possible way to 
discriminate between the two f(R) models is to devise a method able to measure g(z) at 
different scales and redshifts which could, however, be a formidable observational task. On 
the other hand, on the theoretical side, we should correct our estimates of 7 for log k > —1.5 
to take into account deviations from the linear regime 3 which can take place at such large k. 

An alternative parametrization for the growth index has been recently proposed in [87] 
to allow for a redshift dependent growth index. Indeed, one still relies on Eq.(5.4), but now 
7 is assumed to scale with z as 7(2) = 70 + 7a^/(l + z). In order to allow for a larger 
flexibility, we generalize the [87] proposal to 7(2) = 70 + 7az(l + z)~ a and fit the parameters 
(70, 7 a , a) to the g{z) vs £Im(z) relation for each of the point along the chains for both the 
HS and St models. We then plot the median value and the 68% confidence ranges in Fig. 3 
as function of k considering only scales log A; < —1.5 in order to avoid nonlinear effects. As 
a first general comment, we find that the approximation works quite well: it reproduces the 
input g{z) with an error smaller than 1% for all the set of parameters. Moreover, the value 
of a is larger than unity and indicates that the redshift evolution is actually more rapid than 
what is implemented in the [87] parametrization. Although the error bars are significantly 
smaller for the St than for the HS model, the median values of both 70 and j a are very 
similar in the two models up to log k ~ —1.5, while more negative values are obtained for the 
HS model on smaller scales. It is worth stressing that, although small, j a takes a non null 
value and, in particular, \j a \ > 0.02 in agreement with [88] who, based on an analysis of the 
growth of perturbations in GR dark energy models, have argued that |7 a | > 0.02 could be a 
signature for modified gravity. Finally, whereas a ~ 1.735 independent of the scale is a good 

3 One can roughly consider the linear regime valid up to scale where the variance of the matter power 
spectrum is of order unity. For a ACDM model, such an upper limits turns out to be k ~ 0.15/i Mpc -1 , i.e. 
logfc ~ —1.0 for the h ~ 0.7 as we find here. In f(R) gravity, this limit can be also smaller depending on 
which recipe is used to correct the linear power spectrum for the nonlinear effects (which can boost the growth 
of perturbations by several percents) so that, as a conservative choice, we have warn the reader to take with 
caution the results for logfc > —1.5 and exclude them from successive fits. 



-17- 



0.558 ■ 
0.556 
0.554 
=0.552 
0.550 
0.548 
0.546 
0.544 - 



-2.5 -2.0 -1.5 

logk (Mpc -1 ) 



-2.0 -1.5 
log k (Mpc -1 ) 




-2.5 -2.0 -1.5 

log k (Mpc -1 ) 



0.558 ■ 
0.556 
0.554 
0.552 
0.550 
0.548 ■ 



-2.5 -2.0 -1.5 
logk (Mpc -1 ) 



-0.018 r 

-0.020 

-0.022 

,-0.024 

-0.026 

-0.028 

-0.030 



-2.5 -2.0 -1.5 
logk (Mpc -1 ) 



L.90r 
1.85 ; 
1 .80 ; 
I 1.75 : 

i.7o ; 

1.65 : 



-2.5 -2.0 -1.5 
logk (Mpc -1 ) 



Figure 3. Growth index fitting parameters vs logfc for the HS (top) and St (bottom) models. 

approximation for the St model, this is not the case for the HS f(R) theory that shows a clear 
increasing trend with k. Note that, since 7 a ~ —0.028 in both cases, this result show that 
the growth index decreases with z in a faster way for the HS than the St model. Although 
this could suggest a possible way to discriminate between the two models, observationally 
detecting the redshift and scale dependence of the growth index is a rather unrealistic task. 

5.2 The growth factor 

The growth index 7 is not directly observable, but it must be inferred from measurements of 
both the matter density parameter at a given redshift and the growth factor g{z) at the same 
z. This latter quantity may be estimated from redshift space distortions in the galaxy power 
spectra at different z [89, 90]. Actually, what is directly measured is the ratio g{z)/b{z) with 
b(z) the linear bias needed to convert the observed galaxy power spectrum to the matter one. 
If, as a first rough approximation, we assume that the bias is the same for both f{R) and 
ACDM models, we can then solve Eq.(5.3) and compare the predicted growth factor with 
the observed one. 

In order to take care of the uncertainties on the model parameters, we solve Eq.(5.3) 
for all the points in the chain and finally plot the median result having checked that the 95% 
confidence range for g(z) at any given z is negligibly small compared to the statistical uncer- 
tainties. The left panel in Fig. 4 shows the result for the HS model setting k = 0.01 Mpc -1 , 
but the plot is the same as long as we remain in the linear regime (roughly, k < 0.1 Mpc" 1 ). 
We overplot the g(z) data as reported in [91] assembling the different measurements then 
available and here complemented by the recent measurements from the WiggleZ survey [92]. 
As it is apparent, there is a very good match between the model and the data up to z = 1.5, 
while the point at z = 3 (derived from Lya forest power spectrum) appears to disagree. 
However, this latter measurement might be biased by some error in the b(z) estimate at such 
large redshift; in fact, even ACDM, which is successful in many other resepcts, is unable to 
fit this point. Although such an agreement with the growth factor data is surely encouraging, 
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Figure 4. Growth factor as function of z. Left. g(z) vs z for the HS model with k = 0.01 Mpc -1 
with superimposed the observed data. Right. Ratio between the growth factor for the St and HS 
models evaluated at k = 0.01 Mpc" 1 . In both panels, we plot g(z) as inferred from the median of 
the values along the corresponding chains. Note that the spike in the right panel at low z is only a 
numerical artifact and not a physical feature. 



it is worth stressing that they are based on the assumption of same linear bias for both GR 
based ACDM model and modified gravity f{R) scenarios. Therefore, any conclusion drawn 
from Fig. 4 must be taken cum grano salis. 

Finally, we note that the growth factor for the St model also is consistent with the g(z) 
data given the nearly coincidence with the HS model prediction. Indeed, as the right panel 
in Fig. 4 show, the two growth factors closely track each other along the full redshift range 
with a maximum deviation of ~ 1% at very low z which is hard to observationally detect 
with both present and future data. We therefore argue that the growth factor alone can not 
be used to discriminate between the two classes of f(R) theories we are here investigating. 

5.3 The deviations from the GR growth of perturbations 

So far, we have only considered the growth of structures as probed by observable quantities 
related to the matter power spectrum. Actually, the growth of structures can also be probed 
from the point of view of the metric potentials entering the perturbed line element : 

g^dx^dx" = -(1 + 2&)dt 2 + (1 + 2^>)dx 2 . 
While <5 enters the Poisson equation, which in Fourier space reads : 

$k(a) = -4:TTg eff (k,a)(k 2 /a 2 )p m 5 k (a) , 

the other metric potential \& enters the weak lensing potential T = — (<]? — \£)/2 and is related 
to <& through the parameter rj = — ($ + ^)/<I>. In GR, one has Geff/G^ = 1 and r] = (i.e., 
no anisotropic stress), while, for f(R) theories, G e ff(k,a) is given by Eq.(5.2) and 

_ 2(k 2 /a 2 )[r(R)/f(R)] 
V[k ' a) - l + 2{k 2 /a?)[f"{R)/f'{R)\ ■ ^ 

It is then interesting to estimate the theoretical predictions for Q e ff(k, a) and rj(k, a) in order 
to see whether they can help to discriminate between the proposed modified gravity scenarios 
and the concordance ACDM. G e ff an d r/ are shown in Fig. 5 where the central value refers 
to the median and the plotted bar denote 68% confidence range as inferred from the set of 
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Figure 5. Deviations from the GR growth of perturbations in terms of Geff/GN (left) and log?y 
(right) for the HS (top) and St (bottom) models setting k = 0.001 Mpc -1 . 



parameters along the chain. Note that we have set k = 0.001 Mpc , but the results can 
be easily scaled to other values. We also stress that the strong asymmetry of the confidence 
range around the median is a consequence of how Q e ff is defined. Indeed, in the high-z 
regime, both f'(R) and f"(R) become negligibly small and one recovers Geff/GN = 1, while, 
in the low-z limit, all the relevant terms are positive thus leading to Qeff/G^ > 1- As a 
final result, Qeff I G^ can not take values smaller than 1 along the chains thus motivating 
the strongly asymmetric confidence ranges. 

The median values for QeffjG are close to 1 as expected. Indeed, over most of the 
parameter space allowed by the data, both the HS and St Lagrangians are quite similar 
to GR + A so that we indeed expect to recover an effective gravitational constant equal to 
the Newtonian G. Moreover, since the set of parameters have to combine in such a way to 
reproduce the data as well as the ACDM model, the inferred confidence ranges for Q e ff/G 
are quite narrow and deviations from unity become more and more unlikely as z increases 
since, in this regime, both f(R) models reduce to an effective ACDM. A similar qualitative 
discussion also explains the behavior of rj(z) so that deviations of the lensing potential from 
GR are smaller and smaller as z increases. We finally note that, for the St model, Q e ff(k,a) 
and rj(k, a) are always closer to the GR values than for the HS model. We can therefore 
anticipate that the cosmic shear power spectrum is unable to discriminate between the St 
and ACDM models. This is, indeed, what we find in [78] when tomography is not used; on 
the contrary the HS model can give rise to a detectable signature. 
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6 Conclusions 



As soon as the observational and conceptual problems related to the cosmological constant 
and other dark energy scenarios became pressing, modification of the gravity sector of Ein- 
stein field equations immediately appeared as an interesting alternative explanation of the 
observed cosmic accelerated expansion. Fourth -order gravity theories then stand out as the 
most immediate generalization of Einstein GR since they just encoded all the deviations into 
a single analytic function f(R). As the problem of acceleration was solved in this framework, 
a new problem came out, namely how to choose a functional expression for f(R) which is 
able to speed up the expansion and, at the same time, does not violate the local tests of 
gravity and turns off its effect in the early Universe where GR appears to work correctly. 

From a mathematical point of view, one has to look for an f(R) expression satisfying 
some constraints imposed to both fulfill the Solar System tests and recover GR during the 
nucleosynthesis era which is indeed what the HS and St models efficiently do, postulating 
that f(R) is given by Eqs.(2.12) and (2.14), respectively. Our aim here was then to test 
whether these two well behaved models actually fit the data which suggest the accelerated 
expansion of the Universe. To this end, we have considered the background evolution as 
traced by the SNela and GRBs Hubble diagrams, the H(z) data from the differential age 
method applied to passively evolving red galaxies, the BAO constraints on r s {zd)/Dy (z) and 
the WMAP7 distance priors. The combined dataset allows us to probe both the late Universe 
(z < 1.5), the matter dominated era (through GRBs) and the last scattering surface epoch. 
We indeed find that both the HS and St models are able to fit this extended dataset with 
values of (Qm, h, qo,jo) in good agreement with previous results in the literature. We can 
therefore conclude that both the HS and St models are not only theoretically viable f(R) 
proposals, but also observationally motivated choices for a fourth -order gravity Lagrangian. 

An unpleasant outcome of our analysis is represented by the very weak constraints we 
give on the f(R) parameters and the impossibility to discriminate between the two models. 
This is a consequence of the two f(R) functional expressions having been tailored to reduce 
to the ACDM Lagrangian in the high curvature regime. Therefore, the background evolution 
is almost the same and only weak constraints can be set even using extraordinary precise 
data if these only probe the background evolution. Major improvements in both narrowing 
down the confidence ranges and in the possibility to discriminate not only between the HS 
and St models, but, more generally, between dark energy and f(R) theories, are expected 
from the analysis of the growth of perturbations. Indeed, even if one can tailor the f(R) 
parameters in order to closely mimic the same expansion history of a given dark energy 
model, the evolution of the density perturbations can be rather different. In particular, this 
has a strong impact on both the power spectrum and halo statistics [34, 93-96] and the weak 
lensing signals [76-79] and offers the possibility to compare predictions with data and to 
severely constrain the viability of f(R) theories. 

As a preliminary investigation, we have here derived the growth index 7 for the HS 
and St f(R) theories showing that the usual parametrization of the growth rate still holds 
provided a scale dependent 7 is used. However, such a scale dependence is quite weak in 
the deep linear regime (say, e.g., k < 0.01), while the stronger variation for larger k has to 
be confirmed in the nonlinear regime. Moreover, we have also improved the usual growth 
index parametrization allowing for a redshift dependence which turns out to be stronger for 
the HS than for the St model. Unfortunately, using this signature to discriminate between 
the two classes of f(R) theories is theoretically possible, but practically quite difficult (if not 
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impossible) to implement. As an alternative way, one can directly resort to the growth factor 
data as estimated from the redshift distortions of the galaxy power spectrum. Assuming 
that the linear bias is the same as in the GR framework, we have shown that both f(R) 
models agree quite well with the present day data so that it is not possible to discriminate 
among them and GR. Actually, one should first check that the bias is indeed the same in the 
two different frameworks considering that, in f(R) theories, not only the growth of structure 
is different, but also the gravitational potential which plays a key role in the formation of 
galaxies and hence in the determination of their bias function. 

All these preliminary investigations may be considered as a first step towards an im- 
proved analysis of these two classes of f(R) theories. It is indeed the combination of extended 
background evolution tracers (such as SNela and GRBs Hubble diagram, BAO and CMB 
distance priors) and structure growth probes (including galaxy power spectrum and cosmic 
shear) that will finally tell us whether the observed cosmic speed up has been the first evi- 
dence of a new fluid, as mysterious as fascinating, or of new physics in the gravity sector, as 
unexpected as challenging. 
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